HydroNetwork.f90 Source File

!! Manage cells connection in the river network
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!    
!### History
!
! current version  1.0 - 4th October 2016  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      |  04/Oct/2016 | Original code |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Module to manage how cells are connected in a network
!
MODULE HydroNetwork      

! Modules used: 

USE DataTypeSizes, ONLY : &  
! Imported Parameters:
float, short

USE GridLib, ONLY : &
! imported definition:
grid_integer

USE GridOperations, ONLY : &
!Importe routines:
  GetIJ

USE Utilities, ONLY : &
!imported routines:
  GetUnit

USE Loglib, ONLY : &
!Imported routines:
Catch

USE ErrorCodes, ONLY : &
!Imported parameters:
  openFileError

USE StringManipulation, ONLY : &
!Imported routines:
  ToString

IMPLICIT NONE 

! Global (i.e. public) Declarations: 


!TYPE ReachSoilBalance
!	INTEGER :: id
!	!====================================================
!	REAL(KIND = float) :: x0	!!\__beginning of reach
!	REAL(KIND = float) :: y0	!!/                     spatial coordinate
!	REAL(KIND = float) :: x1	!!\__end of reach
!	REAL(KIND = float) :: y1	!!/
!	!====================================================
!	INTEGER :: i0	!!\__beginning of reach
!	INTEGER :: j0	!!/                                 local reference system
!	INTEGER :: i1	!!\__end of reach
!	INTEGER :: j1	!!/
!	!====================================================
!	INTEGER :: ncells !!number of cells in a reach
!	INTEGER :: order !! Horton-Sthraler order
!	REAL(KIND = float) :: slope	!! average reach slope [m/m]
!	REAL(KIND = float) :: length	!!reach length [m]
!	REAL(KIND = float) :: area	!! area drained by end section [m2]
!	INTEGER(KIND = short) :: soilBalanceID
!  REAL(KIND = float) :: n	!!manning roughness coefficient [s m^(-1/3)]
!END TYPE ReachSoilBalance

!TYPE ReachSubsurfaceFlow
!	INTEGER :: id
!	!====================================================
!	REAL(KIND = float) :: x0	!!\__beginning of reach
!	REAL(KIND = float) :: y0	!!/                     spatial coordinate
!	REAL(KIND = float) :: x1	!!\__end of reach
!	REAL(KIND = float) :: y1	!!/
!	!====================================================
!	INTEGER :: i0	!!\__beginning of reach
!	INTEGER :: j0	!!/                                 local reference system
!	INTEGER :: i1	!!\__end of reach
!	INTEGER :: j1	!!/
!	!====================================================
!	INTEGER :: ncells !!number of cells in a reach
!	INTEGER :: order !! Horton-Sthraler order
!	REAL(KIND = float) :: slope	!! average reach slope [m/m]
!	REAL(KIND = float) :: length	!!reach length [m]
!	REAL(KIND = float) :: area	!! area drained by end section [m2]
!	INTEGER            :: routing_model
!END TYPE ReachSubsurfaceFlow

TYPE ReachSurfaceFlow
	INTEGER :: id
	!====================================================
	REAL(KIND = float) :: x0	!!\__beginning of reach
	REAL(KIND = float) :: y0	!!/                     spatial coordinate
	REAL(KIND = float) :: x1	!!\__end of reach
	REAL(KIND = float) :: y1	!!/
	!====================================================
	INTEGER :: i0	!!\__beginning of reach
	INTEGER :: j0	!!/                                 local reference system
	INTEGER :: i1	!!\__end of reach
	INTEGER :: j1	!!/
	!====================================================
	INTEGER :: ncells !!number of cells in a reach
	INTEGER :: order !! Horton-Sthraler order
	REAL(KIND = float) :: slope	!! average reach slope [m/m]
	REAL(KIND = float) :: length	!!reach length [m]
	REAL(KIND = float) :: area	!! area drained by end section [m2]
	INTEGER            :: routing_model
	REAL(KIND = float) :: n	!!manning roughness coefficient [s m^(-1/3)]
	REAL(KIND = float) :: B0	!!bottom width, = 0 for  triangular  section [m]
  REAL(KIND = float) :: alpha !!angle formed by dykes over a horizontal plane [deg]
	REAL(KIND = float) :: k !!flood wave travel time used to compute C1, C2, and C3 [s]
	REAL(KIND = float) :: x !!Muskingum weighting factor [-]
END TYPE ReachSurfaceFlow

TYPE ReachGroundwater
	INTEGER :: id
	!====================================================
	REAL(KIND = float) :: x0	!!\__beginning of reach
	REAL(KIND = float) :: y0	!!/                     spatial coordinate
	REAL(KIND = float) :: x1	!!\__end of reach
	REAL(KIND = float) :: y1	!!/
	!====================================================
	INTEGER :: i0	!!\__beginning of reach
	INTEGER :: j0	!!/                                 local reference system
	INTEGER :: i1	!!\__end of reach
	INTEGER :: j1	!!/
	!====================================================
	INTEGER :: ncells !!number of cells in a reach
	INTEGER :: order !! Horton-Sthraler order
	REAL(KIND = float) :: riverbed ![m]
	REAL(KIND = float) :: scalefactor_conductivity  
END TYPE ReachGroundwater


TYPE ReachSediment
  INTEGER :: id
	!====================================================
	REAL(KIND = float) :: x0	!!\__beginning of reach
	REAL(KIND = float) :: y0	!!/                     spatial coordinate
	REAL(KIND = float) :: x1	!!\__end of reach
	REAL(KIND = float) :: y1	!!/
	!====================================================
	INTEGER :: i0	!!\__beginning of reach
	INTEGER :: j0	!!/                                 local reference system
	INTEGER :: i1	!!\__end of reach
	INTEGER :: j1	!!/
	!====================================================
	INTEGER :: ncells !!number of cells in a reach
	INTEGER :: order !! Horton-Sthraler order
	REAL(KIND = float) :: slope	!! average reach slope [m/m]
	REAL(KIND = float) :: length	!!reach length [m]
	REAL(KIND = float) :: area	!! area drained by end section [m2]
  REAL(KIND = float) :: d50  !!mean sediment size [mm]	
  INTEGER(KIND = short) :: routingMethod
END TYPE ReachSediment

!TYPE NetworkSoilBalance
!	TYPE(ReachSoilBalance),POINTER :: branch (:)
!	INTEGER :: nreach
!END TYPE NetworkSoilBalance

!TYPE NetworkSubsurfaceFlow
!	TYPE(ReachSubsurfaceFlow),POINTER:: branch (:)
!	INTEGER :: nreach
!END TYPE NetworkSubsurfaceFlow

TYPE NetworkGroundwater
	TYPE(ReachGroundwater),POINTER :: branch (:)
	INTEGER :: nreach
END TYPE NetworkGroundwater

TYPE NetworkSurfaceFlow
	TYPE(ReachSurfaceFlow),POINTER :: branch (:)
	INTEGER :: nreach
END TYPE NetworkSurfaceFlow

TYPE NetworkSediment
	TYPE(ReachSediment),POINTER :: branch (:)
	INTEGER :: nreach
END TYPE NetworkSediment


!Public routines
PUBLIC :: ReadHydroNetwork


!Local declarations:


!Local routines
PRIVATE :: ReadNetworkSurfaceFlow
!PRIVATE :: ReadNetworkSubsurfaceFlow
!PRIVATE :: ReadNetworkSoilBalance
PRIVATE :: ReadNetworkGroundwater
PRIVATE :: ReadNetworkSediment

!Local parameters

! Operator definitions:
! Define new operators or overload existing ones.

INTERFACE ReadHydroNetwork
   MODULE PROCEDURE ReadNetworkSurfaceFlow
   !MODULE PROCEDURE ReadNetworkSubsurfaceFlow
   !MODULE PROCEDURE ReadNetworkSoilBalance
   MODULE PROCEDURE ReadNetworkGroundwater
   MODULE PROCEDURE ReadNetworkSediment
END INTERFACE



!=======
CONTAINS
!=======
! Define procedures contained in this module. 
 
  
!==============================================================================
!| Description:
!   Read network for surface flow routing from file
SUBROUTINE ReadNetworkSurfaceFlow &
!
( filename, domain, network )

!arguments with intent in:
CHARACTER (LEN = *), INTENT(IN) :: filename
TYPE(grid_integer), INTENT(IN)  :: domain

!arguments with intent out:
TYPE(NetworkSurfaceFlow), INTENT(OUT) :: network


!local declarations
INTEGER  :: k
INTEGER  :: fileunit
INTEGER  :: err_io
INTEGER  :: id
LOGICAL  :: check_domain

!-------------------------------end of declaration----------------------------- 
!open file in readonly mode
fileunit = GetUnit ()
OPEN (unit = fileunit, file = filename, status = "old", &
     action = "read", iostat = err_io )

IF (err_io /= 0) THEN
  CALL Catch ('error', 'HydroNetwork',        &
              'error opening file of surface flow hydro network: ',    &
              code = openFileError, argument = fileName )
END IF

!count number of branches in  network
k = 0
READ(fileunit,*)
READ(fileunit,*)

DO !loop till end of file	
  READ(fileunit,*,iostat=err_io) id
		IF (err_io /= 0) THEN
			EXIT
    ELSE
      k = k + 1
    END IF
END DO

IF (k < 1) THEN
    CALL Catch ('error', 'HydroNetwork', &
               'no branches detected in surface flow network') 
ELSE
   CALL Catch ('info', 'HydroNetwork', &
               ' number of detected branches in surface flow network: ', &
               argument = ToString(k) ) 
END IF

!rewind file
REWIND (fileunit)
!skip first two lines
READ(fileunit,*)
READ(fileunit,*)

!allocate branches
network % nreach = k
ALLOCATE (network % branch(k))

!read info from file
DO k = 1, network % nreach
  READ (fileunit,*) network % branch (k) % id, network % branch (k) % x0, &
                    network % branch (k) % y0, network % branch (k) % x1, &
                    network % branch (k) % y1, network % branch (k) % ncells, &
                    network % branch (k) % order, network % branch (k) % slope, &
                    network % branch (k) % length, network % branch (k) % area, &
                    network % branch (k) % routing_model, network % branch (k) % n, &
                    network % branch (k) % B0, network % branch (k) % alpha, &
                    network % branch (k) % k, network % branch (k) % x
  !compute i0, j0, i1, j1
  CALL GetIJ (X = network % branch (k) % x0, Y = network % branch (k) % y0, &
              grid = domain, i = network % branch (k) % i0, j = network % branch (k) % j0, &
              check = check_domain)
  IF (.NOT. check_domain) THEN
      CALL Catch ('error', 'HydroNetwork', &
                   'beginning cell out of domain in branch of surface flow with id: ', &
                   argument = ToString(network % branch(k) % id ) )
  END IF

  CALL GetIJ (X = network % branch (k) % x1, Y = network % branch (k) % y1, &
              grid = domain, i = network % branch (k) % i1, j = network % branch (k) % j1, &
              check = check_domain )
  IF (.NOT. check_domain) THEN
      CALL Catch ('error', 'HydroNetwork', &
                   'end cell out of domain in branch of surface flow with id: ', &
                   argument = ToString(network % branch(k) % id ) )
  END IF
END DO

!close file
CLOSE (fileunit)

RETURN
END SUBROUTINE ReadNetworkSurfaceFlow




!==============================================================================
!! Description:
!!   Read network for sub-surface flow routing from file
!SUBROUTINE ReadNetworkSubsurfaceFlow &
!!
!( filename, domain, network )
!
!!arguments with intent in:
!CHARACTER (LEN = *), INTENT(IN) :: filename
!TYPE(grid_integer), INTENT(IN)  :: domain
!
!!arguments with intent out:
!TYPE(NetworkSubsurfaceFlow), INTENT(OUT) :: network
!
!
!!local declarations
!INTEGER  :: k
!INTEGER  :: fileunit
!INTEGER  :: err_io
!INTEGER  :: id
!LOGICAL  :: check_domain
!
!!-------------------------------end of declaration----------------------------- 
!!open file in readonly mode
!fileunit = GetUnit ()
!OPEN (unit = fileunit, file = filename, status = "old", &
!     action = "read", iostat = err_io )
!
!IF (err_io /= 0) THEN
!  CALL Catch ('error', 'HydroNetwork',        &
!              'error opening file of sub-surface flow hydro network: ',    &
!              code = openFileError, argument = filename )
!END IF
!
!!count number of branches in  network
!k = 0
!READ(fileunit,*)
!READ(fileunit,*)
!
!DO !loop till end of file	
!  READ(fileunit,*,iostat=err_io) id
!		IF (err_io /= 0) THEN
!			EXIT
!    ELSE
!      k = k + 1
!    END IF
!END DO
!
!IF (k < 1) THEN
!    CALL Catch ('error', 'HydroNetwork', &
!               'no branches detected in sub-surface flow network') 
!ELSE
!   CALL Catch ('info', 'HydroNetwork', &
!               ' number of detected branches in sub-surface flow network: ', &
!               argument = ToString(k) ) 
!END IF
!
!!rewind file
!REWIND (fileunit)
!!skip first two lines
!READ(fileunit,*)
!READ(fileunit,*)
!
!!allocate branches
!network % nreach = k
!ALLOCATE (network % branch(k))
!
!!read info from file
!DO k = 1, network % nreach
!  READ (fileunit,*) network % branch (k) % id, network % branch (k) % x0, &
!                    network % branch (k) % y0, network % branch (k) % x1, &
!                    network % branch (k) % y1, network % branch (k) % ncells, &
!                    network % branch (k) % order, network % branch (k) % slope, &
!                    network % branch (k) % length, network % branch (k) % area, &
!                    network % branch (k) % routing_model
!  
!  !compute i0, j0, i1, j1
!  CALL GetIJ (X = network % branch (k) % x0, Y = network % branch (k) % y0, &
!              grid = domain, i = network % branch (k) % i0, j = network % branch (k) % j0, &
!              check = check_domain)
!  IF (.NOT. check_domain) THEN
!      CALL Catch ('error', 'HydroNetwork', &
!                   'beginning cell out of domain in branch of sub-surface flow with id: ', &
!                   argument = ToString(network % branch(k) % id ) )
!  END IF
!
!  CALL GetIJ (X = network % branch (k) % x1, Y = network % branch (k) % y1, &
!              grid = domain, i = network % branch (k) % i1, j = network % branch (k) % j1, &
!              check = check_domain )
!  IF (.NOT. check_domain) THEN
!      CALL Catch ('error', 'HydroNetwork', &
!                   'end cell out of domain in branch of sub-surface flow with id: ', &
!                   argument = ToString(network % branch(k) % id ) )
!  END IF
!END DO
!
!!close file
!CLOSE (fileunit)
!
!RETURN
!END SUBROUTINE ReadNetworkSubsurfaceFlow



!==============================================================================
!! Description:
!!   Read network for soil balance from file
!SUBROUTINE ReadNetworkSoilBalance &
!!
!( filename, domain, network )
!
!!arguments with intent in:
!CHARACTER (LEN = *), INTENT(IN) :: filename
!TYPE(grid_integer), INTENT(IN)  :: domain
!
!!arguments with intent out:
!TYPE(NetworkSoilBalance), INTENT(OUT) :: network
!
!
!!local declarations
!INTEGER  :: k
!INTEGER  :: fileunit
!INTEGER  :: err_io
!INTEGER  :: id
!LOGICAL  :: check_domain
!
!!-------------------------------end of declaration----------------------------- 
!!open file in readonly mode
!fileunit = GetUnit ()
!OPEN (unit = fileunit, file = filename, status = "old", &
!     action = "read", iostat = err_io )
!
!IF (err_io /= 0) THEN
!  CALL Catch ('error', 'HydroNetwork',        &
!              'error opening file of soil balance hydro network: ',    &
!              code = openFileError, argument = filename )
!END IF
!
!!count number of branches in  network
!k = 0
!READ(fileunit,*)
!READ(fileunit,*)
!
!DO !loop till end of file	
!  READ(fileunit,*,iostat=err_io) id
!		IF (err_io /= 0) THEN
!			EXIT
!    ELSE
!      k = k + 1
!    END IF
!END DO
!
!IF (k < 1) THEN
!    CALL Catch ('error', 'HydroNetwork', &
!               'no branches detected in soil balance network') 
!ELSE
!   CALL Catch ('info', 'HydroNetwork', &
!               ' number of detected branches in soil balance network: ', &
!               argument = ToString(k) ) 
!END IF
!
!!rewind file
!REWIND (fileunit)
!!skip first two lines
!READ(fileunit,*)
!READ(fileunit,*)
!
!!allocate branches
!network % nreach = k
!ALLOCATE (network % branch(k))
!
!!read info from file
!DO k = 1, network % nreach
!  READ (fileunit,*) network % branch (k) % id, network % branch (k) % x0, &
!                    network % branch (k) % y0, network % branch (k) % x1, &
!                    network % branch (k) % y1, network % branch (k) % ncells, &
!                    network % branch (k) % order, network % branch (k) % slope, &
!                    network % branch (k) % length, network % branch (k) % area, &
!                    network % branch (k) % soilBalanceID, network % branch (k) % n
!  
!  !compute i0, j0, i1, j1
!  CALL GetIJ (X = network % branch (k) % x0, Y = network % branch (k) % y0, &
!              grid = domain, i = network % branch (k) % i0, j = network % branch (k) % j0, &
!              check = check_domain)
!  IF (.NOT. check_domain) THEN
!      CALL Catch ('error', 'HydroNetwork', &
!                   'beginning cell out of domain in branch of soil balance with id: ', &
!                   argument = ToString(network % branch(k) % id ) )
!  END IF
!
!  CALL GetIJ (X = network % branch (k) % x1, Y = network % branch (k) % y1, &
!              grid = domain, i = network % branch (k) % i1, j = network % branch (k) % j1, &
!              check = check_domain )
!  IF (.NOT. check_domain) THEN
!      CALL Catch ('error', 'HydroNetwork', &
!                   'end cell out of domain in branch of soil balance with id: ', &
!                   argument = ToString(network % branch(k) % id ) )
!  END IF
!END DO
!
!!close file
!CLOSE (fileunit)
!
!RETURN
!END SUBROUTINE ReadNetworkSoilBalance



!==============================================================================
!| Description:
!   Read network for groundwater-surface interaction from file
SUBROUTINE ReadNetworkGroundwater &
!
( filename, domain, network )

!arguments with intent in:
CHARACTER (LEN = *), INTENT(IN) :: filename
TYPE(grid_integer), INTENT(IN)  :: domain

!arguments with intent out:
TYPE(NetworkGroundwater), INTENT(OUT) :: network


!local declarations
INTEGER  :: k
INTEGER  :: fileunit
INTEGER  :: err_io
INTEGER  :: id
LOGICAL  :: check_domain

!-------------------------------end of declaration----------------------------- 
!open file in readonly mode
fileunit = GetUnit ()
OPEN (unit = fileunit, file = filename, status = "old", &
     action = "read", iostat = err_io )

IF (err_io /= 0) THEN
  CALL Catch ('error', 'HydroNetwork',        &
              'error opening file of groundwater hydro network: ',    &
              code = openFileError, argument = filename )
END IF

!count number of branches in  network
k = 0
READ(fileunit,*)
READ(fileunit,*)

DO !loop till end of file	
  READ(fileunit,*,iostat=err_io) id
		IF (err_io /= 0) THEN
			EXIT
    ELSE
      k = k + 1
    END IF
END DO

IF (k < 1) THEN
    CALL Catch ('error', 'HydroNetwork', &
               'no branches detected in groundwater network') 
ELSE
   CALL Catch ('info', 'HydroNetwork', &
               ' number of detected branches in groundwater network: ', &
               argument = ToString(k) ) 
END IF

!rewind file
REWIND (fileunit)
!skip first two lines
READ(fileunit,*)
READ(fileunit,*)

!allocate branches
network % nreach = k
ALLOCATE (network % branch(k))

!read info from file
DO k = 1, network % nreach
  READ (fileunit,*) network % branch (k) % id, network % branch (k) % x0, &
                    network % branch (k) % y0, network % branch (k) % x1, &
                    network % branch (k) % y1, network % branch (k) % ncells, &
                    network % branch (k) % order, network % branch (k) % riverbed, &
                     network % branch (k) % scalefactor_conductivity 
  !compute i0, j0, i1, j1
  CALL GetIJ (X = network % branch (k) % x0, Y = network % branch (k) % y0, &
              grid = domain, i = network % branch (k) % i0, j = network % branch (k) % j0, &
              check = check_domain)
  IF (.NOT. check_domain) THEN
      CALL Catch ('error', 'HydroNetwork', &
                   'beginning cell out of domain in branch of groundwater with id: ', &
                   argument = ToString(network % branch(k) % id ) )
  END IF

  CALL GetIJ (X = network % branch (k) % x1, Y = network % branch (k) % y1, &
              grid = domain, i = network % branch (k) % i1, j = network % branch (k) % j1, &
              check = check_domain )
  IF (.NOT. check_domain) THEN
      CALL Catch ('error', 'HydroNetwork', &
                   'end cell out of domain in branch of groundwater with id: ', &
                   argument = ToString(network % branch(k) % id ) )
  END IF
END DO

!close file
CLOSE (fileunit)

RETURN
END SUBROUTINE ReadNetworkGroundwater


!==============================================================================
!| Description:
!   Read network for sediment routing from file
SUBROUTINE ReadNetworkSediment &
!
( filename, domain, network )

!arguments with intent in:
CHARACTER (LEN = *), INTENT(IN) :: filename
TYPE(grid_integer), INTENT(IN)  :: domain

!arguments with intent out:
TYPE(NetworkSediment), INTENT(OUT) :: network


!local declarations
INTEGER  :: k
INTEGER  :: fileunit
INTEGER  :: err_io
INTEGER  :: id
LOGICAL  :: check_domain

!-------------------------------end of declaration----------------------------- 
!open file in readonly mode
fileunit = GetUnit ()
OPEN (unit = fileunit, file = filename, status = "old", &
     action = "read", iostat = err_io )

IF (err_io /= 0) THEN
  CALL Catch ('error', 'HydroNetwork',        &
              'error opening file of sediment routing hydro network: ',    &
              code = openFileError, argument = filename )
END IF

!count number of branches in  network
k = 0
READ(fileunit,*)
READ(fileunit,*)

DO !loop till end of file	
  READ(fileunit,*,iostat=err_io) id
		IF (err_io /= 0) THEN
			EXIT
    ELSE
      k = k + 1
    END IF
END DO

IF (k < 1) THEN
    CALL Catch ('error', 'HydroNetwork', &
               'no branches detected in sediment routing network') 
ELSE
   CALL Catch ('info', 'HydroNetwork', &
               ' number of detected branches in sediment routing network: ', &
               argument = ToString(k) ) 
END IF

!rewind file
REWIND (fileunit)
!skip first two lines
READ(fileunit,*)
READ(fileunit,*)

!allocate branches
network % nreach = k
ALLOCATE (network % branch(k))

!read info from file
DO k = 1, network % nreach
  READ (fileunit,*) network % branch (k) % id, network % branch (k) % x0, &
                    network % branch (k) % y0, network % branch (k) % x1, &
                    network % branch (k) % y1, network % branch (k) % ncells, &
                    network % branch (k) % order, network % branch (k) % slope, &
                    network % branch (k) % length, network % branch (k) % area, &
                    network % branch (k) % d50
  
  !compute i0, j0, i1, j1
  CALL GetIJ (X = network % branch (k) % x0, Y = network % branch (k) % y0, &
              grid = domain, i = network % branch (k) % i0, j = network % branch (k) % j0, &
              check = check_domain)
  IF (.NOT. check_domain) THEN
      CALL Catch ('error', 'HydroNetwork', &
                   'beginning cell out of domain in branch of sediment transport with id: ', &
                   argument = ToString(network % branch(k) % id ) )
  END IF

  CALL GetIJ (X = network % branch (k) % x1, Y = network % branch (k) % y1, &
              grid = domain, i = network % branch (k) % i1, j = network % branch (k) % j1, &
              check = check_domain )
  IF (.NOT. check_domain) THEN
      CALL Catch ('error', 'HydroNetwork', &
                   'end cell out of domain in branch of sediment transport with id: ', &
                   argument = ToString(network % branch(k) % id ) )
  END IF
END DO

!close file
CLOSE (fileunit)

RETURN
END SUBROUTINE ReadNetworkSediment




END MODULE HydroNetwork